\documentclass[12pt, a4paper]{article}
\usepackage{geometry}
\geometry{
	a4paper,
	left=12.7 mm,
	right=12.7 mm,
	top=12.7 mm,
	bottom=12.7 mm,
}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{physics}
\graphicspath{ {./} }
\usepackage{ctex}

\newcommand{\bvec}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\formula}[1]{\text{式} \ref{#1} }
\begin{document}

\section{插值问题}
\footnote{本文是Sauer《数值分析》，欧阳杰等人《数值分析》的学习笔记}插值问题要解决一个很简单、也很现实的问题：
我们想找出一个函数 $y=f(x)$， 经过（或者说，靠近\footnote{部分教材可能限定插值函数必须经过所有给定的点}）一组已知的$N$个点：
\begin{equation} \label{eq_1}
\{(x_1, y_1),(x_2, y_2),...,(x_N, y_N)\}
\end{equation}
以下，我们假定插值函数 $y=f(x)$是一个具有有限项的多项式函数：
\begin{equation}  \label{eq_2}
	f(x) = a_0 + a_1 x + a_2 x^2+...+a_n x^n
\end{equation}
问题转换为，找到这个多项式的所有系数$\{a_0, a_1, a_2,..., a_n\}$。

$f(x)$的阶数（最高次项的次数$x^n$）取决于具体的方法。
不过可以肯定的是，\textbf{$N$个点最多插值出$N-1$阶的多项式；
	并且，如果插值得到了$N-1$阶的多项式，那他应该经过所有给定的点}。
如果试图用$N$个点插值$N$阶及以上多项式，会发现相应的高阶项系数为$0$。

\subsection{最小二乘法}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.4\linewidth]{interpolate}
	\caption{随机生成4个点，并分别进行$1,2,3$阶插值。注意到$3$阶插值函数经过了全部$4$个点}
	\label{fig:interpolate}
\end{figure}


结合 \formula{eq_1} 与 \formula{eq_2}，我们得到一个线性方程组：
\begin{equation}
	\left \{
	\begin{aligned}
		a_0 + a_1 x_1 + a_2 x_1^2 + ... &= y_0\\
		a_0 + a_1 x_2 + a_2 x_2^2 + ... &= y_1\\
		a_0 + a_1 x_3 + a_2 x_3^2 + ... &= y_2\\
		&...
	\end{aligned}
	\right.
\end{equation}
此处$x_0, x_1,...$是已知的、反而$a_0, a_1,...$是未知的。
问题变成求解这个线性方程组的问题。
将其改写为矩阵形式：
\begin{equation}
	\left (
	\begin{matrix}
		1 & x_1 & x_1^2 & ... \\
		1 & x_2 & x_2^2 & ... \\
		1 & x_3 & x_3^2 & ... \\
		... & 
	\end{matrix}
	\right)
	\left (
	\begin{matrix}
		a_0&\\
		a_1&\\
		a_2&\\
		... &\\
	\end{matrix}
	\right)
	=
	\left (
	\begin{matrix}
		y_1&\\
		y_2&\\
		y_3&\\
		... &\\
	\end{matrix}
	\right)
\end{equation}
或
\begin{equation}
	X a = y
\end{equation}
求解$a$，好像只需
$$
a = X^{-1} y
$$
不过，我们有点太理想了：$X$很可能是不可逆的。不过没关系，我们有\textbf{最小二乘法}：
\begin{equation}
X^T X a = X^T y \Rightarrow a = (X^T X)^{-1} X^T y
\end{equation}
至于为什么可以这么做就说来话长了，\textsl{读者自证不难}。总而言之，最小二乘法是一种既简单又有效的方法。
如果数据量不大，即使只使用最基本的Octave/Matlab代码，也能在几行之内搞定。
事实上，Octave/Matlab提供了polyfit和polyval函数来简化多项式拟合。

对于最小二乘法，$f(x)$的阶数是我们自定的，从一阶的线性函数到$N-1$阶的多项式函数均可。这给了我们很大的灵活性。

\subsection{Lagrange插值法}
Lagrange插值法可以用于$N-1$阶的插值。Lagrange法的插值结果是
\begin{equation}
	f(x) = \sum_i^N l_i(x) y_i
\end{equation}
其中$l_i(x)$是Lagrange基函数
\begin{equation}
	l_i(x) = \frac{\Pi_{j\ne i} (x-x_j)}{\Pi_{j\ne i} (x_i-x_j)} = \frac{(x-x_1)(x-x_2)...(x-x_{i-1})(x-x_{i+1})...(x-x_N)}{(x_i-x_0)(x_i-x_1)...(x_i-x_{i-1})(x_i-x_{i+1})...(x_i-x_N)}
\end{equation}
因为$j=1,2,...,N$，但是没有$j=i$的项，因此$l_i$是$N-1$阶的、$f(x)$也是$N-1$阶的

\subsection{Newton插值法（插商法）}
Newton插值法也是$N-1$阶的插值。Newton法的插值结果是
\begin{equation}
	\begin{aligned}
		f(x) &= \sum_{i=1}^N f[x_1, x_2, ..., x_i] \Pi_{j=1}^{i-1} (x-x_j)\\
		& = f[x_1] + f[x_1, x_2] (x-x_1) + f[x_1, x_2, x_3] (x-x_1) (x-x_2) +...\\
		& + f[x_1, x_2, ..., x_N]  (x-x_1) (x-x_2) ... (x-x_{N-1}) 
	\end{aligned}
\end{equation}
其中 $f[x_1, x_2, ..., x_k] $ 称为插商，定义如下：
\begin{itemize}
	\item $f[x_1] = y_1$
	\item $f[x_1, x_2] = \frac{y_2 - y_1}{x_2 - x_1}$
	\item $f[x_1, x_2, x_3] = \frac{f[x_2, x_3]-f[x_1, x_2]}{x_3 - x_1}$
	\item $f[x_1, x_2, x_3, x_4] = \frac{f[x_2, x_3, x_4]-f[x_1, x_2, x_3]}{x_4 - x_1}$
	\item ...
\end{itemize}
可以画出插商表以帮助计算插商。Newton法的优势是，很容易新增一个点。

如果我们分别使用最小二乘法、Newton法、Lagrange法得到$N-1$阶的多项式，那这些多项式有什么区别吗？有哪一个更准确的问题吗？
答案是\textbf{没有。这三种方法插值得到的$N-1$阶多项式是一样的，并且经过所有给定的点。}这无疑是一个好消息。

\section{分段插值}
有时我们会遇到插值具有成百上千个点的问题。
技术上说，使用一个成百上千阶的多项式是可行的，
但这么高阶的多项式往往导致严重的\textbf{过拟合}，以及引入大量意义不明的拟合参数。
每一本数值分析书都会涉及这个问题，因此按下不表。

这个问题无法被\textsl{解决}，但可以考虑一些替代方案：
要么使用最小二乘法、使用所有点拟合一个低阶多项式；
要么将点分为数个小组、在每一个小组上分别拟合低阶多项式。
比如说，在 $(x_1, x_3), (x_4, x_6), (x_7, x_9),...$等数个小区间上分别拟合$2$阶多项式。


\section{带导数的插值问题}

现在，我们不仅知道一组点，甚至还知道一组导数。共有$M$个这样的条件：
\begin{equation} \label{eq_3}
	\{(x_1, y_1),(x_2, y_2),...,(x_N, y_N),(x_{N+1}, y'_{N+1}),(x_{N+2}, y'_{N+2}), ...,(x_M, y'_M)\}
\end{equation}
在这种情况下，我们也同样可以构造最多$M-1$阶的多项式。

\subsection{最小二乘法}
我仍然推荐最小二乘法：
对\formula{eq_2}逐项求导：
\begin{equation}  \label{eq_4}
	f'(x) = a_1 x + 2 a_2 x+ 3 a_3 x^2 + ...+n a_n x^{n-1}
\end{equation}
我们发现，导数的条件相当于引入了额外的一组线性方程；
结合 \formula{eq_2}, \formula{eq_3},  \formula{eq_4}，现在的线性方程组是：
\begin{equation}
	\left (
	\begin{matrix}
		1 & x_1 & x_1^2 &... \\
		1 & x_2 & x_2^2 &... \\
		1 & x_3 & x_3^2 &... \\
		... &  \\
		0 & 1 & 2x_{N+1} &... \\
		0 & 1 & 2x_{N+2} &... \\
		... & \\
	\end{matrix}
	\right)
	\left (
	\begin{matrix}
		a_0&\\
		a_1&\\
		a_2&\\
		... &\\
	\end{matrix}
	\right)
	=
	\left (
	\begin{matrix}
		y_1&\\
		y_2&\\
		y_3&\\
		... &\\
		y'_{N+1}  & \\
		y'_{N+2} & \\
		... &\\
	\end{matrix}
	\right)
\end{equation}
仍然可以用最小二乘法的思路解决。

\subsection{拼凑法}
有时，导数的条件比较少，并且刚好同时知道某处的函数值和导数值。比如，我们假定知道$x_1$处的导数。这样，我们一共有$M=N+1$个条件：
\begin{equation} 
	\{(x_1, y_1),(x_2, y_2),...,(x_{N}, y_{N}), (x_1, y'_1)\}
\end{equation}
并且，对于$\{(x_1, y_1),(x_2, y_2),...,(x_{N}, y_{N})\}$这$N$个点，我们已经插值得到了一个$N-1$阶的多项式 $f_{N-1}(x)$。
为了引入导数条件，我们实则没有必要从头来过。我们设最终的$M-1=N+1-1=N$阶插值函数是
\footnote{请注意这里的符号问题。不同文章因为符号约定不同会有点差异。}
\begin{equation}
	\begin{aligned}
		f(x) &= f_{N-1}(x) + c \Pi_{j = 1}^{N} (x-x_j) \\
		& = f_{N-1}(x) + c (x-x_1)(x-x_2)...(x-x_{N}) \\
	\end{aligned}
\end{equation}
对其求导并代入 $x=x_1$（这里有一个可以大幅简化计算的小技巧）：
\begin{equation}
	\begin{aligned}
		f'(x_1) = f'_{N-1}(x_1) + c (x_1-x_2)...(x_1-x_{N}) \\
	\end{aligned}
\end{equation}
由于$f'(x_1) = y'_1$是已知的，我们由此确定$c$
$$
c = \frac{y'_1 - f'_{N-1}(x_1)}{(x_1-x_2)...(x_1-x_{N})}
$$
以及最终的$f(x)$。

\section{自定义基的插值}
\textsl{jojo，我不用多项式啦}。我们也可以使用一组自定义的函数（比如，三角函数）
\begin{equation}
	\{f_0(x), f_1(x), f_2(x),..., f_n(x)\}
\end{equation}
充当插值函数 $y=f(x)$的基：
\begin{equation} 
	f(x) = a_0  f_0(x) + a_1 f_1(x)  + a_2 f_2(x) +...+a_n f_n(x) 
\end{equation}
这种情况下，我们仍然可以列出线性方程组：
\begin{equation}
	\left (
	\begin{matrix}
		f_0(x_1) & f_1(x_1) & f_2(x_1)  & ... \\
		f_0(x_2)  & f_1(x_2)  & f_2(x_2) & ... \\
		f_0(x_3)  & f_1(x_3) & f_2(x_3)  & ... \\
		... & 
	\end{matrix}
	\right)
	\left (
	\begin{matrix}
		a_0&\\
		a_1&\\
		a_2&\\
		... &\\
	\end{matrix}
	\right)
	=
	\left (
	\begin{matrix}
		y_1&\\
		y_2&\\
		y_3&\\
		... &\\
	\end{matrix}
	\right)
\end{equation}
并且沿用最小二乘法的思路解决它。某种意义上，没有太多区别。
此时前一个矩阵也称为Gram矩阵。

\end{document}
